##########
# R Info #
##########
version
#platform       x86_64-w64-mingw32          
#arch           x86_64                      
#os             mingw32                     
#system         x86_64, mingw32             
#status                                     
#major          3                           
#minor          6.0                         
#year           2019                        
#month          04                          
#day            26                          
#svn rev        76424                       
#language       R                           
#version.string R version 3.6.0 (2019-04-26)
#nickname       Planting of a Tree 

library(ggplot2)
library(ggpubr)
library(maps)
library(wesanderson)
library(stargazer)
library(dplyr)

sessionInfo()
#R version 3.6.0 (2019-04-26)
#Platform: x86_64-w64-mingw32/x64 (64-bit)
#Running under: Windows 10 x64 (build 19041)
#Matrix products: default
#locale:
#[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
#[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
#[5] LC_TIME=English_United States.1252    
#attached base packages:
#[1] stats     graphics  grDevices utils     datasets  methods   base     
#other attached packages:
#[1] dplyr_1.0.4       stargazer_5.2.2   wesanderson_0.3.6 maps_3.3.0        ggpubr_0.4.0     
#[6] ggplot2_3.3.2    
#loaded via a namespace (and not attached):
#[1] zip_2.0.4         Rcpp_1.0.3        cellranger_1.1.0  pillar_1.4.3      compiler_3.6.0   
#[6] forcats_0.5.0     tools_3.6.0       lifecycle_0.2.0   tibble_3.0.4      gtable_0.3.0     
#[11] pkgconfig_2.0.3   rlang_0.4.10      openxlsx_4.1.4    DBI_1.1.0         rstudioapi_0.11  
#[16] curl_4.3          haven_2.2.0       rio_0.5.16        rJava_0.9-11      withr_2.4.1      
#[21] dplyr_1.0.4       hms_0.5.3         generics_0.0.2    xlsxjars_0.6.1    vctrs_0.3.6      
#[26] grid_3.6.0        tidyselect_1.1.0  data.table_1.12.8 glue_1.4.1        R6_2.4.1         
#[31] rstatix_0.6.0     readxl_1.3.1      foreign_0.8-71    carData_3.0-3     car_3.0-7        
#[36] purrr_0.3.3       tidyr_1.1.2       magrittr_1.5      backports_1.1.5   scales_1.1.0     
#[41] ellipsis_0.3.0    abind_1.4-5       assertthat_0.2.1  colorspace_1.4-1  ggsignif_0.6.0   
#[46] xlsx_0.6.1        stringi_1.4.5     munsell_0.5.0     broom_0.7.5       crayon_1.3.4




#############
# Load Data #
#############
load("PostDataJQD.RData")

colnames(df)
# country:          country name
# year:             year of posting
# yearmonth:        year + month of posting
# party_country:    party name + country name
# party:            party name
# cmp_familyid:     party family id in the CMP
# cmp_family:       party family name in the CMP
# ches_ideology:    general left-right ideology in the CHES
# gps_populism:     populism score in the GPS
# gps_econlr:       economic left-right ideology in the GPS      
# gps_soclr:        social left-right ideology in the GPS
# Page.Name:        name of facebook page
# Likes.at.Posting: number of page likes at posting
# Likes:            number of Likes
# Comments:         number of comments
# Shares:           number of shares
# Love:             number of Love reactions
# Wow:              number of Wow reactions
# Haha:             number of Haha reactions 
# Sad:              number of Sad reactions
# Angry:            number of Angry reactions
# totalreactions:   sum of Likes, Love, Wow, Haha, Sad, and Angry
# loveprop:         proportion of Love reactions (Love / totalreactions)
# angryprop:        proportion of Angry reactions (Angry / totalreactions)

summary(df)

# number of countries
length(unique(df$country)) # 79

# number of parties
length(unique(df$party_country)) # 690

# compute emotional polarization
df$lapolarization <- (df$angryprop + 1) / (df$loveprop + 1)

# set colors
rgy <- wes_palette("Darjeeling1", 3)




##############
# Appendix C #
##############
# like proportion
df$likeprop <- df$Likes / df$totalreactions
df$likeprop[df$totalreactions == 0] <- 0  

# country level data
countries <- data.frame(country=sort(unique(df$country)))

# average number of total reactions and like proportion by country
countries$totalreactionmean <- tapply(df$totalreactions, df$country, mean)
countries$avelikeprop <- tapply(df$likeprop, df$country, mean)

# function to get confidence intervals
getCI <- function(vector){
  m <- mean(vector)
  s <- sd(vector)
  error <- qt(0.995,df=length(vector)-1) * s / sqrt(length(vector))
  left <- m - error
  right <- m + error
  return(c(left, right))
}

#pdf("ReactionLikeRank.pdf", width=12, height=15)
par(mfrow=c(1,2), mar=c(3.1,8,3.1,2.1))

# order countries by average total reactions
countries <- countries[order(countries$totalreactionmean),]

# figure C1 left panel
plot(countries$totalreactionmean, 1:79, type="n",
     ylab="", yaxt="n", ylim=c(3,77),
     xlab="", xlim=c(0,2400),
     main="Average # Reactions per Post")
abline(h=1:79, col="gray95")
abline(v=mean(df$totalreaction), lty=2)
points(countries$totalreactionmean, 1:79, pch=20, cex=2)
axis(2, 1:79, labels=countries$country, las=1, cex.axis=0.8)
text(mean(df$totalreaction),1, labels="Global Mean", cex=0.7)
for(i in 1:79){
  ci <- getCI(df$totalreaction[df$country == countries$country[i]])
  segments(ci[1], i, ci[2], i, lwd=1.3)
}

# order countries by like proportion
countries <- countries[order(countries$avelikeprop),]

# figure C1 right panel
plot(countries$avelikeprop, 1:79, type="n",
     ylab="", yaxt="n", ylim=c(3,77),
     xlab="", xlim=c(0.7,1),
     main="Average Like Proportion")
abline(h=1:79, col="gray95")
abline(v=mean(df$likeprop), lty=2)
points(countries$avelikeprop, 1:79, pch=20, cex=2)
axis(2, 1:79, labels=countries$country, las=1, cex.axis=0.8)
text(mean(df$likeprop),1, labels="Global Mean", cex=0.7)
for(i in 1:79){
  ci <- getCI(df$likeprop[df$country == countries$country[i]])
  segments(ci[1], i, ci[2], i, lwd=1.3)
}

#dev.off()


# calculate reaction diversity by post
df$wowprop <- df$Wow / df$totalreactions
df$wowprop[df$totalreactions == 0] <- 0

df$hahaprop <- df$Haha / df$totalreactions
df$hahaprop[df$totalreactions == 0] <- 0

df$sadprop <- df$Sad / df$totalreactions
df$sadprop[df$totalreactions == 0] <- 0

df$reactionfrac <- apply(df[,c("loveprop", "angryprop", "likeprop",
                               "wowprop", "hahaprop", "sadprop")], 1, function(p){
                                 1 - sum(p^2)
                               })
df$reactionfrac[which(df$totalreactions == 0)] <- NA
summary(df$reactionfrac)

# group posts by total reaction quantile
fracq <- data.frame(group=1:10,
                    quantile=seq(0.1, 1, by=0.1))

df$totalreactions_log <- log(df$totalreactions)

df$fracgroup <- NA
for(i in 10:1){
  rn <- which(df$totalreactions_log <= quantile(df$totalreactions_log, 
                                                fracq$quantile[i],
                                                na.rm=TRUE))
  df$fracgroup[rn] <- fracq$group[i]
}

# figure C2
#pdf("fraction.pdf", width=7, height=5)
par(mfrow=c(1,1), mar=c(4.1,4.1,2.1,2.1))

boxplot(df$reactionfrac ~ df$fracgroup,
        outline=FALSE,
        ylab="Reaction Fractionalization",
        xlab="Quantile of Total Number of Reactions")

#dev.off()




##############
# Appendix D #
##############
# country level data
countries <- data.frame(country=sort(unique(df$country)))

# compute correlation between love and angry proportions
countries$loveangrycorr <- sapply(countries$country, function(p){
  cor(df$loveprop[df$country == p], df$angryprop[df$country == p])
})

# order countries by correlation
countries <- countries[order(countries$loveangrycorr),]

# figure D1
#pdf("CorrelationLoveAngry.pdf", width=7, height=15)
par(mfrow=c(1,1), mar=c(3.1,8,3.1,2.1))

plot(countries$loveangrycorr, 1:79, type="n",
     ylab="", yaxt="n", ylim=c(3,77),
     xlab="", xlim=c(min(countries$loveangrycorr)-0.03, max(countries$loveangrycorr)+0.03),
     main="Correlation between Love and Angry Proportions")
abline(h=1:79, col="gray95")
abline(v=0, lty=2)
points(countries$loveangrycorr, 1:79, pch=18, cex=1.8)
axis(2, 1:79, labels=countries$country, las=1, cex.axis=0.8)

#dev.off()




##############
# Appendix E #
##############
threedatesets <- unique(df[,c("party_country", "country",
                              "cmp_family", "ches_ideology", 
                              "gps_populism", "gps_econlr", "gps_soclr")])

# table E1
table(threedatesets$cmp_family)


# table E2
stargazer(threedatesets[,c("ches_ideology", 
                           "gps_econlr", "gps_soclr", "gps_populism")],
          covariate.labels=c("CHES General Ideology", 
                             "GPS Economic Ideology", "GPS Social Ideology", 
                             "GPS Populism"),
          omit.summary.stat=c("p25","p75"),  digits=2)


# correlations
round(cor(threedatesets[,c("ches_ideology", 
                           "gps_econlr", "gps_soclr", "gps_populism")],
          use="pairwise.complete.obs"), 2)


# figure E1
#pdf("partydistribution.pdf", width=9, height=5.5)
par(mfrow=c(2,2))

plot(density(threedatesets$ches_ideology, na.rm=TRUE),
     main="(1) CHES General Ideology")
rug(threedatesets$ches_ideology, col=rgb(0.1,0.1,0.1,0.5))

plot(density(threedatesets$gps_econlr, na.rm=TRUE),
     main="(2) GPS Economic Ideology")
rug(threedatesets$gps_econlr, col=rgb(0.1,0.1,0.1,0.5))

plot(density(threedatesets$gps_soclr, na.rm=TRUE),
     main="(3) GPS Social Ideology")
rug(threedatesets$gps_soclr, col=rgb(0.1,0.1,0.1,0.5))

plot(density(threedatesets$gps_populism, na.rm=TRUE),
     main="(4) GPS Populism")
rug(threedatesets$gps_populism, col=rgb(0.1,0.1,0.1,0.5))

#dev.off()




##############
# Appendix F #
##############
# set memory limit
memory.limit(500000000)

# calculate Love proportion, Angry proportion, and emotional polarization
# using Likes.at.Posting as a denominator
df$loveprop2 <- (df$Love / df$Likes.at.Posting) * 100
summary(df$loveprop2)

df$angryprop2 <- (df$Angry / df$Likes.at.Posting) * 100
summary(df$angryprop2)

df$lapolarization2 <- (df$angryprop2 + 1) / (df$loveprop2 + 1)
summary(df$lapolarization2)

# set ggplot theme
plain2 <- theme(
  panel.border = element_rect(fill = NA),
  panel.grid =  element_line(colour = "grey92"),
  strip.background = element_rect(fill = "grey85", 
                                  colour = "grey20"),
  panel.background = element_rect(fill = "white"),
  plot.title = element_text(hjust = 0, face = "bold"),
  legend.position = "none",
  axis.title.x = element_text(size = 15),
  axis.title.y = element_text(size = 15)
)


# party family
partyfamily <- unique(df[,c("cmp_familyid", "cmp_family")])
partyfamily <- partyfamily[-1,]
partyfamily <- partyfamily[order(partyfamily$cmp_familyid),]

partyfamily$`Party Family` <- factor(partyfamily$cmp_family, 
                                     levels=partyfamily$cmp_family)

partyfamily <- rbind(partyfamily, partyfamily)

partyfamily$Reaction <- rep(c("Love", "Angry"), each=10)
partyfamily$Reaction <- factor(partyfamily$Reaction, levels=c("Love", "Angry"))

# calculate average love and angry proportions by party family
partyfamily$prop[1:10] <- tapply(df$loveprop2, df$cmp_familyid, mean, na.rm=T)
partyfamily$prop[11:20] <- tapply(df$angryprop2, df$cmp_familyid, mean, na.rm=T)

partyfamily$xloc <- 1:10

partyfamily$type <- "Average Love and Angry Proportions"

partyfamily2 <- cbind(partyfamily[1:10,1:5], partyfamily[11:20,5])
colnames(partyfamily2)[5:6] <- c("loveprop", "angryprop")

partyfamily2$xloc <- 1:10

# figure F1 left panel
party_diff <- ggplot(data = partyfamily, 
                     aes(x = xloc, y = prop, colour = Reaction)) + 
  facet_wrap(. ~ type, ncol = 1) +
  geom_segment(data = partyfamily2, aes(x = xloc, xend = xloc, y = loveprop, yend = angryprop),
               lwd=0.8, color = "gray60") +
  geom_point(cex=4, pch=19, alpha=0.7) +
  scale_color_manual(breaks = c("Love", "Angry"), 
                     values = rgy[2:1]) +
  scale_x_discrete(limits=factor(1:10), labels= partyfamily$cmp_family[1:10]) + 
  xlab("") +
  ylab("Love and Angry Proportion") +
  coord_cartesian(ylim=c(0, 0.0004*100)) +
  theme_bw() +
  theme(axis.text.x=element_text(color = "black", size=11, angle=30, vjust=.8, hjust=0.8),
        #panel.grid = element_blank(),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        legend.title = element_text(size =13),
        legend.text = element_text(size =13),
        legend.position = c(0.12, 0.12),
        legend.background = element_rect(size=0.7, linetype="solid", colour="black"))


# calculate average emotional polarization by party
partyfamily2$lapolarization <- (partyfamily2$angryprop + 1) / (partyfamily2$loveprop + 1)

partyfamily2 <- partyfamily2[order(partyfamily2$lapolarization, decreasing=FALSE),]
partyfamily2$xloc <- 1:10

partyfamily2$type <- "Emotional Polarization"

# figure F1 right panel
party_ep <- ggplot(data = partyfamily2, 
                   aes(x = xloc, y = lapolarization)) + 
  geom_hline(yintercept=1, lty=2) +
  geom_line(color = "gray60", lwd=0.8) +
  geom_point(cex=6, pch=18, colour=rgy[3]) +
  facet_wrap(. ~ type, ncol = 1) +
  scale_x_discrete(limits=factor(1:10), labels= partyfamily2$cmp_family) + 
  xlab("") +
  ylab("Emotional Polarization") +
  theme_bw() +
  theme(axis.text.x=element_text(color = "black", size=11, angle=30, vjust=.8, hjust=0.8),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        legend.title = element_text(size =13),
        legend.text = element_text(size =13))


# figure F1 all
ggarrange(party_diff, party_ep,
          ncol = 2, nrow = 1)
#ggsave("PartyFamilyPageLike.pdf", width=12, height=6)


# CHES general ideology
set.seed(123456)
psamples <- sample(intersect(which(!is.na(df$ches_ideology)),
                             which(!is.na(df$Likes.at.Posting))), 20000)

il <- df[psamples, c("ches_ideology", "loveprop2")]
colnames(il)[2] <- 'prop'
il$type <- "Love"

ia <- df[psamples, c("ches_ideology", "angryprop2")]
colnames(ia)[2] <- "prop"
ia$type <- "Angry"

ip <- df[psamples, c("ches_ideology", "lapolarization2")]
colnames(ip)[2] <- "prop"
ip$type <- "Emotional Polarization"

df2 <- rbind(il, ia, ip)

df2$type <- factor(df2$type, levels=c("Love", "Angry", "Emotional Polarization"))

tapply(df2$prop, df2$type, min)
tapply(df2$prop, df2$type, max)

blank_data <- data.frame(type=rep(unique(df2$type), each=2), 
                         ches_ideology=0, 
                         prop=c(0,0.00065*100,0,0.00065*100,0.98,1.03))
df2$rugloc <- rep(c(0,0,0.98), each=nrow(df2)/3)

ches <- ggplot() +
  geom_blank(data = blank_data, aes(x = ches_ideology, y = prop)) +
  geom_hline(data = data.frame(yint=1, type="Emotional Polarization"), 
             aes(yintercept = yint), lty=2) +
  geom_smooth(data = df2, 
              aes(x = ches_ideology, y = prop, colour = type),
              method = "loess", lwd = 1.2,
              level=0.99) +
  facet_wrap(. ~ type, ncol = 3, scales = "free") +
  scale_color_manual(breaks = c("Love", "Angry", "Emotional Polarization"), 
                     values = rgy[c(2:1,3)]) +
  ggtitle(label = "A. CHES General Ideology") +
  xlab("Left---Right") + 
  ylab("") +
  plain2 + 
  scale_y_continuous(breaks = seq(0,2,by=0.0003*100))


# GPS economic ideology
set.seed(123456)
psamples <- sample(intersect(which(!is.na(df$gps_econlr)),
                             which(!is.na(df$Likes.at.Posting))), 20000)

el <- df[psamples, c("gps_econlr", "loveprop2")]
colnames(el)[2] <- 'prop'
el$type <- "Love"

ea <- df[psamples, c("gps_econlr", "angryprop2")]
colnames(ea)[2] <- "prop"
ea$type <- "Angry"

ep <- df[psamples, c("gps_econlr", "lapolarization2")]
colnames(ep)[2] <- "prop"
ep$type <- "Emotional Polarization"

df2 <- rbind(el, ea, ep)

df2$type <- factor(df2$type, levels=c("Love", "Angry", "Emotional Polarization"))

blank_data <- data.frame(type=rep(unique(df2$type), each=2), 
                         gps_econlr=0, 
                         prop=c(0,0.0011*100,0,0.0011*100,0.93,1.03))
df2$rugloc <- rep(c(0,0,0.93), each=nrow(df2)/3)

econ <- ggplot() +
  geom_blank(data = blank_data, aes(x = gps_econlr, y = prop)) +
  geom_hline(data = data.frame(yint=1, type="Emotional Polarization"), 
             aes(yintercept = yint), lty=2) +
  geom_smooth(data = df2, 
              aes(x = gps_econlr, y = prop, colour = type),
              method = "loess", lwd = 1.2, 
              level=0.99) +
  facet_wrap(. ~ type, ncol = 3, scales = "free") +
  scale_color_manual(breaks = c("Love", "Angry", "Emotional Polarization"), 
                     values = rgy[c(2:1,3)]) +
  ggtitle(label = "B. GPS Economic Ideology") +
  xlab("Pro State---Pro Market") + 
  ylab("") +
  plain2 + 
  scale_y_continuous(breaks = seq(0,2,by=0.0003*100))


# GPS social ideology
set.seed(123456)
psamples <- sample(intersect(which(!is.na(df$gps_soclr)), 
                             which(!is.na(df$Likes.at.Posting))), 20000)

sl <- df[psamples, c("gps_soclr", "loveprop2")]
colnames(sl)[2] <- 'prop'
sl$type <- "Love"

sa <- df[psamples, c("gps_soclr", "angryprop2")]
colnames(sa)[2] <- "prop"
sa$type <- "Angry"

sp <- df[psamples, c("gps_soclr", "lapolarization2")]
colnames(sp)[2] <- "prop"
sp$type <- "Emotional Polarization"

df2 <- rbind(sl, sa, sp)

df2$type <- factor(df2$type, levels=c("Love", "Angry", "Emotional Polarization"))

blank_data <- data.frame(type=rep(unique(df2$type), each=2), 
                         gps_soclr=0, 
                         prop=c(0,0.00035*100,0,0.00035*100,0.99,1.02))
df2$rugloc <- rep(c(0,0,0.99), each=nrow(df2)/3)

soc <- ggplot() +
  geom_blank(data = blank_data, aes(x = gps_soclr, y = prop)) +
  geom_hline(data = data.frame(yint=1, type="Emotional Polarization"), 
             aes(yintercept = yint), lty=2) +
  geom_smooth(data = df2, 
              aes(x = gps_soclr, y = prop, colour = type),
              method = "loess", lwd = 1.2, 
              level=0.99) +
  facet_wrap(. ~ type, ncol = 3, scales = "free") +
  scale_color_manual(breaks = c("Love", "Angry", "Emotional Polarization"), 
                     values = rgy[c(2:1,3)]) +
  ggtitle(label = "C. GPS Social Ideology") +
  xlab("Liberal---Conservative") + 
  ylab("") +
  plain2 + 
  scale_y_continuous(breaks = seq(0,2,by=0.0003*100))


# GPS populism
set.seed(123456)
psamples <- sample(intersect(which(!is.na(df$gps_populism)),
                             which(!is.na(df$Likes.at.Posting))), 20000)

pl <- df[psamples, c("gps_populism", "loveprop2")]
colnames(pl)[2] <- 'prop'
pl$type <- "Love"

pa <- df[psamples, c("gps_populism", "angryprop2")]
colnames(pa)[2] <- "prop"
pa$type <- "Angry"

pp <- df[psamples, c("gps_populism", "lapolarization2")]
colnames(pp)[2] <- "prop"
pp$type <- "Emotional Polarization"

df2 <- rbind(pl, pa, pp)

df2$type <- factor(df2$type, levels=c("Love", "Angry", "Emotional Polarization"))

blank_data <- data.frame(type=rep(unique(df2$type), each=2), 
                         gps_populism=0, 
                         prop=c(0,0.00045*100,0,0.00045*100,0.98,1.02))
df2$rugloc <- rep(c(0,0,0.98), each=nrow(df2)/3)

pop <- ggplot() +
  geom_blank(data = blank_data, aes(x = gps_populism, y = prop)) +
  geom_hline(data = data.frame(yint=1, type="Emotional Polarization"), 
             aes(yintercept = yint), lty=2) +
  geom_smooth(data = df2, 
              aes(x = gps_populism, y = prop, colour = type),
              method = "loess", lwd = 1.2,
              level=0.99) +
  facet_wrap(. ~ type, ncol = 3, scales = "free") +
  scale_color_manual(breaks = c("Love", "Angry", "Emotional Polarization"), 
                     values = rgy[c(2:1,3)]) +
  ggtitle(label = "D. GPS Populism") +
  xlab("Populism") + 
  ylab("") +
  plain2 + 
  scale_y_continuous(breaks = seq(0,2,by=0.0003*100))


# figure F2 all
ggarrange(ches, econ, soc, pop,
          ncol = 1, nrow = 4)
#ggsave("IdeologyPageLike.pdf", width=12, height=12)




##############
# Appendix G #
##############
# load country-month and party-month data
load("CountryMonthTrend.RData")
load("PartyMonthTrend.RData")

# set colors
rgy <- wes_palette("Darjeeling1", 3)

# use 10 countries
# 10 countries with many reactions excluding US
countries10 <- c("Belgium", "Canada", "Netherlands", "Australia", "Argentina", 
                 "Poland", "Spain", "Japan", "Czech Republic", "New Zealand")


# create country-month data
cm2 <- country_month[country_month$country %in% countries10,]

colnames(cm2)

cml <- cm2[,c(1,2,3,9)]
colnames(cml)[4] <- "prop"
cml$Reaction <- "Love"

cma <- cm2[,c(1,2,3,10)]
colnames(cma)[4] <- "prop"
cma$Reaction <- "Angry"

cm3 <- rbind(cml, cma)

cm3$countryreaction <- paste0(cm3$country, " (", cm3$Reaction, ")")
cm3$countryreaction <- factor(cm3$countryreaction, 
                              unique(cm3$countryreaction[order(cm3$country)]))

cm3$yearmonthorder <- as.numeric(as.factor(cm3$yearmonth))


# create party-month data
pm2 <- party_month[party_month$country %in% countries10,]

colnames(pm2)

pml <- pm2[,c(1,2,4,9)]
colnames(pml)[4] <- "prop"
pml$Reaction <- "Love"

pma <- pm2[,c(1,2,4,10)]
colnames(pma)[4] <- "prop"
pma$Reaction <- "Angry"

pm3 <- rbind(pml, pma)

pm3$countryreaction <- paste0(pm3$country, " (", pm3$Reaction, ")")
pm3$countryreaction <- factor(pm3$countryreaction, 
                              unique(pm3$countryreaction[order(pm3$country)]))

pm3$yearmonthorder <- as.numeric(as.factor(pm3$yearmonth))


# create election-month data
el <- data.frame(country=c("Argentina", "Australia", "Australia",
                           "Belgium", "Canada", "Czech Republic",
                           "Japan", "Netherlands", "New Zealand",
                           "Poland", "Spain", "Spain", "Spain"),
                 yearmonth=c("2019-10", "2016-07", "2019-05",
                             "2019-05", "2019-10", "2017-10",
                             "2017-10", "2017-03", "2017-09",
                             "2019-10", "2016-06", "2019-04", "2019-11"),
                 stringsAsFactors=FALSE)

e2 <- rbind(el, el)

e2$Reaction <- rep(c("Love", "Angry"), each=nrow(e2)/2)

e2$countryreaction <- paste0(e2$country, " (", e2$Reaction, ")")

e2$countryreaction <- factor(e2$countryreaction, 
                             unique(e2$countryreaction[order(e2$country)]))

e2$yearmonth <- factor(e2$yearmonth, levels=levels(country_month$yearmonth))
e2$yearmonthorder <- as.numeric(e2$yearmonth)


# figure G1
ggplot() + 
  geom_vline(xintercept = c(10.5,22.5,34.5,46.5), lty=1, col = "gray90", lwd = 0.3) +
  geom_vline(data = e2, aes(xintercept = yearmonthorder), lty=2) +
  facet_wrap(. ~ countryreaction, ncol = 2) +
  scale_x_continuous(breaks = c(10.5,22.5,34.5,46.5), labels = 2017:2020) +
  geom_line(data = pm3, 
            aes(x = yearmonthorder, y = prop, group = party), 
            lwd=0.5, colour="gray80") +
  geom_line(data = cm3, 
            aes(x = yearmonthorder, y = prop, colour = Reaction), lwd=1.3) + 
  scale_color_manual(breaks = c("Love", "Angry"), 
                     values = rgy[c(2:1)]) +
  xlab("March 2016-February 2020") + 
  ylab("Love and Angry Proportion") +
  coord_cartesian(ylim=c(0, 0.5)) +
  theme_bw() +
  theme(legend.position = "none",
        panel.grid = element_blank(),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15))
#ggsave("MonthTrend.pdf", width=12, height=15)
